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We derive and analyze in the framework of the mild-slope approximation a new 
double-layer Boussinesq-type model which is linearly and nonlincarly accurate up 
to deep water. Assuming the flow to be irrotational, we formulate the problem 
in terms of the velocity potential thereby lowering the number of unknowns. The 
model derivation combines two approaches, namely the method proposed by Agnon 
et al. (Agnon et al. 1999 J. Fluid Mech. 399, 319-333) and enhanced by Madsen 
et al. (Madsen et al. 2003 Proc. R. Soc. Lond. A 459, 1075-1104) which consists 
in constructing infinite-series Taylor solutions to the Laplace equation, to trun- 
cate them at a finite order and to use Pade approximants, and the double-layer 
approach of Lynett & Liu (Lynett & Liu 2004 Proc. R. Soc. Lond. A 460, 2637- 
2669) allowing to lower the order of derivatives. We formulate the model in terms 
of a static Dirichlet-Neumann operator translated from the free surface to the still- 
water level, and we derive an approximate inverse of this operator that can be built 
once and for all. The final model consists of only four equations both in one and 
two horizontal dimensions, and includes only second-order derivatives, which is a 
major improvement in comparison with so-called high-order Boussincsq models. A 
linear analysis of the model is performed and its properties arc optimized using a 
free parameter determining the position of the interface between the two layers. 
Excellent dispersion and shoaling properties are obtained, allowing the model to 
be applied up to the deep-water value kh = 10. Finally, numerical simulations are 
performed to quantify the nonlinear behaviour of the model, and the results exhibit 
a nonlinear range of validity reaching at least kh = 
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1. Introduction 

During the past two decades, the original Boussinesq (1872) model for flat bottom 
and its uneven bottom version derived by Peregrine (1967) have been widely studied 
and extended to tackle more and more realistic physical problems. Consequently, 
Boussinesq-type models have emerged as an attractive and commonly used tool 
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for coastal applications and engineering purposes. The derivation of such models 
is based on a polynomial approximation of the vertical profile of the velocity field, 
which allows to reduce the problem size from three to two space dimensions. These 
models are usually formulated as conservation equations for mass and momentum, 
including spatial and temporal derivatives of the free surface elevation and the 
velocity. In practice, their range of applicability is measured in terms of kh, where 
k is the wavenumber and h the water depth. 

The conventional Boussinesq model for uneven bottoms (Peregrine 1967), which 
employs a quadratic polynomial approximation for the vertical flow distribution, 
is a depth- averaged model based on two fundamental assumptions, namely weak 
nonlincarity and weak dispersion. Its range of applicability is limited to kh < 0.75, 
as stated in Madsen et al. (2002, 2003), so that this model has poor dispersion 
properties in intermediate depths. Moreover, the weakly nonlinear assumption lim- 
its the largest wave height that can be modelled accurately. As a result, substantial 
efforts have been devoted to extend the linear and nonlinear range of applicability 
of Boussinesq-type models. The first historical improvement, proposed by Nwogu 
(1993), consists in using a reference velocity at a specified depth, allowing the re- 
sulting model to be linearly applicable at intermediate depths. Similar models for 
short-amplitude and long waves have been more recently proposed and rigorously 
justified (Bona et al. 2002, 2005; Chazel 2007). An effort similar to the one of Nwogu 
(1993) was pursued by Madsen & S0rensen (1992), which was followed by the works 
of Liu (1994) and Wei et al. (1995) in which the authors efficiently removed the 
weak nonlincarity assumption, allowing the model to simulate wave propagation 
with strong nonlinear interaction. According to the reviews proposed by Madsen & 
Schaffer (1998, 1999), these new Boussinesq-type models allow to extend the linear 
range of applicability to kh = 6, but it turned out that a similar improvement on 
the nonlinear characteristics was very difficult to reach. Then, so-called high-order 
Boussinesq-type models were derived to further enhance the deep-water linear and 
nonlinear accuracy using a higher-order (at least fourth-order) polynomial approx- 
imation for the vertical flow distribution. One such example is the formulation of 
Gobbi et al. (2000) which uses a fourth-order polynomial approximation: this model 
exhibits excellent linear properties up to kh = 6 for the dispersion relation and up 
to kh = 4 for the vertical profiles of orbital velocities, whereas nonlinear behaviour 
is fairly well captured up to kh = 3. The price to pay for such an improvement is a 
significant increase in computational cost, since this model includes up to fifth-order 
derivatives and therefore requires a complex numerical scheme. 

Over the past decade, two parallel approaches have emerged, one aiming at ex- 
tending even more the range of applicability of the model of Gobbi et al. (2000) into 
very deep water without increasing too much the numerical complexity, and another 
aiming at lowering the numerical cost of the latter while at least preserving the lin- 
ear and nonlinear properties. The first approach has been extensively investigated 
by Madsen and co-workers with a first breakthrough in 1999 (Agnon et al. 1999). 
In this work, the authors presented a new procedure by which it was possible to 
achieve the same accuracy on both linear and nonlinear properties. The main idea is 
to obtain approximate solutions to the Laplace equation (combined with the exact 
nonlinear free surface and bottom conditions) through truncated scries expansions. 
While the formulation of Agnon et al. (1999) involves velocity variables evaluated 
at the still-water level and is limited to kh = 6, the extended method proposed 
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by Madscn et al. (2002, 2003) completely removes the conventional shallow-water 
limitation, allowing for modelling fully nonlinear and highly dispersive waves up 
to kh w 40, i.e. in very deep water. This extended approach is based on velocity 
variables taken at optimized levels and optimal expansions through the use of Pade 
approximants. An extension to rapidly varying bathymetry has been proposed re- 
cently (Madsen et al. 2006). Although a few numerical approaches based on this 
model have been proposed (Fuhrman & Bingham 2004; Engsig-Karup et al. 2006, 
2008), the counterpart for this wide range of applicability is the numerical com- 
plexity of the underlying model which includes derivatives up to fifth-order and 
consists of more equations (and more unknowns) than the alternative model of 
Gobbi et al. (2000). An interesting alternative approach has been chosen by Jamois 
et al. (2006), where the authors propose to use a velocity potential formulation and 
to truncate the infinite series expansions of the solutions to the Laplace equation 
at a lower order: the resulting model is linearly and nonlinearly accurate only up to 
kh = 10, but entails a much lighter numerical complexity with less equations and 
with derivatives up to fourth-order only. 

The second approach that has been studied is the double-layer approach, as pro- 
posed among others by Lynett & Liu (2004a, b) and Audusse (2005). This approach 
is based on the idea of trading fewer unknowns and higher spatial derivatives for 
more unknowns and lower spatial derivatives. The multi-layering concept developed 
in the above references can be seen as an efficient mathematical tool to reduce the 
order of derivatives in any model, while increasing its linear and nonlinear range of 
applicability. The double-layer modelling proposed by Lynett & Liu (2004a, b) is 
purely conceptual since the two layers have the same density. However, the result- 
ing model, which is based on classical depth-integrated Boussinesq-type equations, 
allows to model accurately wave propagation up to kh = 6, both linearly and non- 
linearly. This offers a very interesting alternative to high-order models, since this 
double-layer model is less complex (including derivatives up to third-order only) 
and more accurate in deep water. 

The present work is mainly inspired by these two approaches, namely the one 
of Madsen et al. (2002, 2003) and the one of Lynett & Liu (2004a, 6). The primary 
goal of this paper is to offer an efficient alternative to the models of Madsen et al. 
(2002, 2003) by mixing their procedure, the simplifications of Jamois et al. (2006), 
and the double-layer concept of Lynett & Liu (2004a, b). We aim at deriving a 
model which is 1) applicable to complex domains such as coastal areas, islands 
or estuaries and 2) linearly and nonlinearly accurate up to deep water, but with 
lower complexity than the previous models (i.e. lower order of derivatives and lower 
number of equations). The model proposed herein satisfies all these conditions as it 
exhibits excellent linear and nonlinear dispersive properties up to kh = 10, consists 
of four equations in both one and two horizontal dimensions (denoted by 1DH and 
2DH, respectively), and includes up to second-order spatial derivatives only. 

The present model hinges on a static Dirichlct-Ncumann operator and its ap- 
proximation using a double-layer technique and Pade approximants. The advantage 
of using the static operator (that is, defined on a fixed domain) as opposed to the 
usual Dirichlet-Neumann operator defined at the free surface is that the static 
operator (or its approximation) can be computed once and for all. The Dirichlet- 
Neumann operator has been extensively investigated over the past two decades. 
Exact expressions of the static operator, thereby leading to exact dispersion rela- 
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tions, can be found in the work of Craig & Sulem (1993), Dommermuth & Yue 
(1987), Smith (1998), and in the work of Matsuno (1993) and its extension to very 
general bathymetries by Artiles & Nachbin (2004a, b). However, the application of 
the above approaches to complex 2DH domains has not been reported yet. Thus, 
with an eye toward coastal engineering applications, we prefer to base our approach 
on approximating the static Dirichlet-Neumann operator. Furthermore, we mention 
the new promising approach of Ablowitz et al. (2006) based on an exact integral 
representation of the usual Dirichlet-Neumann operator where no approximation 
is needed. Finally, we observe that the double-layer technique used to derive the 
approximate static Dirichlet-Neumann operator helps reducing the order of the 
derivatives while improving the accuracy of the model. 

The paper is organized as follows. In §2, our model is formulated in terms of 
a static Dirichlet-Neumann operator, and we derive in §3 an approximation to 
this operator. A linear analysis of the model is presented in §4, including linear 
dispersion, vertical profiles of velocities, and linear shoaling. These properties are 
optimized based on Stokes linear wave theory, and it is shown that the model is 
accurate even for deep water conditions. Finally, in §5, numerical simulations are 
developed in 1DH to assess the nonlinear properties of the model for flat bottom 
conditions. 

2. Derivation of the double-layer formulation 

(a) Governing equations and boundary conditions 

We aim at formulating a double-layer Boussincsq-type model for the three- 
dimensional irrotational flow of an inviscid and incompressible fluid with a free 
surface. We focus here on so-called gravity waves or water waves, i.e. the evolution of 
a fluid under the only influence of gravity. The capillary effects due to the presence of 
surface tension arc neglected. Moreover, we assume constant atmospheric pressure 
at the free surface of the fluid. We adopt a Cartesian coordinate system, where we 
denote by X — (x, y) the horizontal coordinates and by z the vertical one, the z-axis 
pointing upwards. The time-dependent fluid domain is bounded from below by the 
static sea bottom and from above by the time-dependent free surface. We restrict 
this study to the case where the bathymetry and the free surface elevation arc 
single-valued continuous functions, i.e. they can be described by the graph of two 
functions X i— > z(X) = —h(X) and (t,X) i— > n(t,X) respectively. The level z = 
corresponds to the still- water level. As shown in figure 1, the fluid is divided into 
two layers by an interface z = z(X) = —ah(X) where a is an arbitrary parameter 
in ]0, 1[. Thus, the thickness of the two layers are constant fractions of the still- 
water depth and do not depend on the free surface elevation. Unless the bottom 
is flat, the interface level z is therefore spatially (but not temporally) variable. 
The upper layer of the fluid is denoted by and the lower layer by f^i namely 
^{t) = { (X, z) ; z(X) < z < r](t, X)} and ft 2 = { {X, z) ; z(X) < z < z{X)}. We 
point out that this fluid division into two layers is purely conceptual since both 
layers have the same density. 

As far as the bathymetry is concerned, we assume in this work that the still- 
water depth h verifies |V/i| <C 1, which corresponds to the classical mild-slope 
approximation. This approximation consists in neglecting all the quadratic (and 

Article submitted to Royal Society 



A double-layer Boussinesq-type model 



■5 




X=(x,y) 

Figure 1. Representation of the fluid domain 



higher order) terms in V/i along with the derivatives of h of order greater than 
one. Physically, this approximation means that we assume the wavelength of the 
free surface waves to be shorter than the distance over which the bathymetry (and 
thus the still- water depth) varies appreciably. We point out that, in the mild-slope 
framework, the overall amplitude of bottom topography levels can still be large. 
The mild-slope approximation plays an important role in the derivation and lin- 
ear optimization of the present model, and it seems quite arduous to incorporate 
higher-order bathymetric terms without significantly increasing the complexity of 
the model. For very general bathymetries in 1DH, we refer to Artiles & Nachbin 
(2004a, b). For clarity, all the equations derived using this mild-slope approximation 
are indicated in this work by the symbol f=a instead of the equality symbol. 

Since the flow is assumed irrotational, there exists a velocity potential </> such 
that u = V</>, if = d z <j>, where u denotes the horizontal velocity of the fluid, 
w the vertical velocity, and V the horizontal gradient operator (d x , d y ) T . We 
define the velocity potentials <f>i and the vertical velocities Wi in each layer by 
<f>i = 0^ , Wi = d z 4>i where the subscript i 6 {1, 2} denotes the layer index. 

The fluid motion in each layer is governed by the following equations written in 
terms of the velocity potential 4>i and the vertical velocity Wi , 

A& + fi£& = 0, (X,z)eSk, (2-la) 

dtfa + i|V0 4 | 2 + X -w\ + gz + Pi =R, (X,z)G fli , (2.1b) 

where Pi is the reduced pressure field in each layer, g the gravitational acceleration, 
and R the Bernoulli constant. Equation (2.1a) corresponds to the Laplace (or con- 
tinuity) equation and (2.1b) corresponds to the Bernoulli (or momentum) equation. 
The Bernoulli constant R only depends on time. Therefore, up to a time-dependent 
shift of the velocity potential, we can take this constant equal to P a tm where P a tm 
is the constant atmospheric pressure at the free surface. 

At the free surface z = r/(t,X), the following boundary conditions are enforced: 

/ d t -q + V-q ■ V0i - ioi = , at z = r) , (2.2a) 

I Pi = P atm , at Z = 7? , (2.2b) 

where (2.2a) is the usual kinematic free surface condition expressing that the free 
surface is a bounding surface, i.e. no fluid particle can cross it. At the interface z = 
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z(X) between the layers, the following natural continuity conditions are enforced: 
4>i = 4>2 , W\=W2, P\= f 2 , at z = z . (2.3) 

Observe that <p\ = 4> 2 and w\ = w 2 at z = z imply V0i = V</>2 and thus u\ = u 2 
at z = z, hence recovering the continuity conditions on the horizontal and vertical 
velocities enforced by Lynett & Liu (2004a). 

Finally, the system of equations is closed by the usual kinematic condition at the 
sea bottom z = z(X): 

V/i • V<j> 2 + w 2 = , &tz = z, (2.4) 
which expresses that the sea bottom is a bounding surface. We now introduce 

' Mt, X) = (j>x{t, X,z = r)(t, X)), w Y = Wl (t, X,z = r)(t, X)) , 
< $ i (t,X) = ct) i {t,X,z = z(X)), w~i=Wi(t,X,z = z{X)) , 

t Mt, X) = X, z = z{X)) , W2- = w 2 {t,X,z = z(X)) . 

Following Zakharov (1968), we can reformulate the equations as 

dtJi + ^|V^| 2 - i^ 2 (l + |V?7| 2 ) + gr, = , (2.5a) 

9^ + V77-V^ 1 -5JI(l + |V77| 2 )=0, (2.5b) 

where (2.5a) is the Euler equation expressed at the free surface and (2.5b) the 
kinematic condition at the free surface, 

A0i + d z wi = , z < z < r\ , (2.5c) 

A0 2 + d z w 2 = , z < z < z , (2.5d) 

where (2.5c) and (2.5d) are the Laplace equations in each layer, and 



where (2.5e) and (2.5f) correspond to the continuity conditions at the interface, 
and (2.5g) to the kinematic condition at the bottom, where we used the mild- 
slope hypothesis to neglect the |V/i| 2 term. We point out that we work with a 
velocity potential formulation, unlike Madsen et al. (2002, 2003) who formulated 
the governing equations in terms of the velocity variables u and w. This choice stems 
from our will to minimize the total number of equations in the model: in 2DH and 
in the present double-layer framework under the irrotational flow assumption, our 
velocity potential formulation allows to consider two less equations than with a 
velocity formulation. 

In the sequel, equations (2.5a) and (2.5b) are left unchanged as they define the 
fully nonlinear time-stepping problem. We focus on the Laplace equations and the 
remaining boundary conditions to close the time-stepping problem by expressing the 



Article submitted to Royal Society 



01 = 4>2 , 
Wl = W2 , 

wi + Vh- VTfo w , 



(2.5e) 
(2.5f) 
(2-5g) 
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vertical velocity at the free surface uT\ in terms of the velocity potential (j>\ at that 
surface, the free surface 77, and the bathymetry h. This relation corresponds to the 
well-known Dirichlct-Neumann operator. The next subsections and §3 are devoted 
to the crucial construction of an accurate, yet computationally cheap approximation 
to this operator. 



(6) A translated Dirichlet- Neumann operator 

The Dirichlet-Neumann operator associated to the problem (2.5c)-(2.5g) is de- 
noted by Q[rj, h] and is defined by Q[rj, h]tp — d z 4>i\ z=ri for any smooth enough 
function ip, where (^1,(^2) solves the boundary value problem composed of equa- 
tions (2.5c)-(2.5g) along with the Dirichlet condition <\>\ = tp at z = rj. One can 
thus simply write the closure between the unknowns w\ , <j>i , and 77 as 

v% = g[ri,h]fa. (2.6) 

This Dirichlct-Neumann operator is at the heart of the derivation of Boussinesq- 
type models since the structure and accuracy of these models essentially depend on 
the method used to construct an approximation to this operator. Once this approx- 
imation is derived, there are two options. The first one is to eliminate the vertical 
velocity variable w{ from the equations by plugging (2.6) into (2.5a), (2.5b). This 
method has been classically used in Boussinesq-type models and has the advantage 
of lowering the number of equations to solve at each time step, but significantly 
increases their complexity. The other option has been used for instance by Madsen 
et al. (2002, 2003) and consists in keeping wi in the equations, which entails to solve 
(2.5a), (2.5b) and then to compute w{ through the use of the Dirichlet-Neumann 
operator Q [77, h] at each time step. This is the method we have chosen to use to 
keep equations complexity to a minimum. 

The main difficulty in finding an approximation to the Dirichlet-Neumann op- 
erator is that it involves solving the Laplace equations (2.5c) and (2.5d) along 
with the boundary conditions (2.5e)-(2.5g) and (f>i = 4>\ at z = 77, on a time- 
dependent domain bounded from above by the free surface z = rj. Keeping w{ in 
the equations involves constructing an approximation to Q[rj, h] at each time step, 
which can be a serious computational problem as it increases the numerical cost. 
An interesting work-around to this issue consists in constructing an alternative 
Dirichlet-Neumann operator expressed at the still-water level, and then finding a 
closure between the unknown functions at the free surface z = 77 and the ones at 
the still- water level z = 0. To this end, we introduce 4>o(t, X) = 4>\{t, X, z = 0) and 
wo(t,X) = wi(t,X,z = 0) and denote by Go[h] the Dirichlct-Neumann operator 
corresponding to 77 = 0, i.e. 

g [h] = g[o,h}. (2.7) 

The translated operator Go[h] is such that go[h]tp = d z 4>i\z=a = where {(t>\,(j>2) 
solves the boundary value problem 

A<t3i + d 2 z (j)i = , z < z < , 

Acj) 2 + dlfo = , z < z < z , 

4>i=ip, atz = 0, (2.8) 

<t>i = <i>2 , d z <j)i = d z <p 2 , at z = z , 

d z (f>2 + V/i • V02 = , at z = z . 
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The main advantage of this translated Dirichlct-Ncumann operator is that the 
boundary value problem (2.8) is posed on a static domain bounded from above by 
the still-water level z = and from below by the sea bottom z = z. 

Remark 2.1. We point out that 4>q is not defined everywhere since the bounding 
free surface r\ can take negative values. This issue can be solved by extending the 
solutions to the Laplace equation above z = rj when n < 0, using the fact that both <f>\ 
and its normal derivative w\ are continuous at the free surface. This mathematical 
trick allows to artificially define 0o andwo when r] < 0. This tool has been implicitly 
used by many authors such as Nwogu (1993), Wei et al. (1995), Gobbi et al. (2000) 
or Madsen et al. (2002, 2003) who derived models based on an horizontal velocity 
(or potential) variable taken at free- surface independent levels, thus allowing these 
levels to exceed the bounding value z — rj for r) negative enough. 



(c) Closure relation and model formulation 

Our goal is twofold. Firstly to construct an approximation to the translated 
Dirichlct-Neumann operator ^o[^]j and secondly to look for closure relations be- 
tween the unknowns <f>i, w{, 4>q, and Wq. The second objective can be readily 
achieved via a Taylor expansion of <pi and W\ at the still-water level z = 0. In- 
deed, combining the MacLaurin expansions of <j)\ and w\ at respectively the fourth 
and third order (see Remark 2.2) and the Laplace equation A<j>i = —d z wi at z = 
yields the desired closure relations, namely 



J 0i = (1 - OuA) <f> + (fii - 71 A) w , 
\wi = -Pi A 0o + (1 - 51 A) w , 

r/ 2 — _ rf 

where a\ = —, pi = r), and 71 = — . 

We can now state our double-layer model as follows: 



(2.9) 



|V0i| 2 



1 



5Ji(l + |V77| 2 ) + . 977 = 



d t n + V77 • V0i - wi{l + IV77I 2 ) = , 



_ vA + {1 _IL A) g«pp [h] 



= 4>i , 



(2.10) 



where Gq PP [h] is an approximation to the static Dirichlet-Neumann operator Qq [h] 
which is detailed in the following section. 



The main advantages of this model are that 1) it consists of only four equa- 
tions both in 1DH and 2DH, 2) it can be used in complex domains, and 3) the 
approximate Dirichlet-Neumann operator Qq VP [h] can be computed once and for 
all at t = since this operator is static. Furthermore, we will see in §3 that it only 
includes at most second-order horizontal derivatives. This is a major improvement 
in comparison with high-order Boussinesq-type models such as those of Jamois et 
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al. (2006) and Madsen et al. (2002, 2003) which contain respectively fourth- and 
fifth-order derivatives, and consist of respectively five equations in 1DH and 2DH, 
and five equations in f DH and seven in 2DH. 

Remark 2.2. The truncation orders used in (2.9) can be motivated by a dimen- 
sional analysis. We scale the vertical coordinate z and the free surface rj with the 
typical amplitude a of the surface waves, the horizontal coordinate X by the typical 
wavelength A and introduce the mean depth h$. The truncation errors in the two 
equations of (2.9) are respectively of order O^ji 2 ,e 5 /i 2 ) and 0(e 4 /x 2 , e 3 // 2 ), where 
e = a I ho and /i = ft. 2 ,/ A 2 correspond respectively to the nonlinearity and dispersion 
parameters. Thus, the truncated terms are almost third and fourth powers of the 
steepness parameter s defined by s = e^/JL, whose maximum value is admittedly 
Smax ~ 0.142 for a stable wave (Williams 1981). Combining this value and the fact 
that we consider fully nonlinear waves, for which e is of order 0(1), motivates the 
truncation order in (2.9). 

3. An approximate static Dirichlet-Neumann operator 

(a) Theoretical solutions to the Laplace equations 

The first step in the derivation of the approximate static Dirichlet-Neumann 
operator Go PP [h] is to look for solutions to the Laplace equations 

Acf>i + dl4>i = , (X,z)€SU, (3.1) 

where we have redefined the upper-layer domain as fii = { (X, z) ; z(X) < z < 0}. 
To this end, we follow the generalized Boussinesq procedure introduced by Madsen 
et al. (2002, 2003) which consists in looking for a solution under the form of an 
infinite Taylor series in the vertical coordinate. The main difference between this 
method and the classical Boussinesq procedure is that in the latter, one looks for 
a finite series solution in the vertical coordinate, i.e. a low-order polynomial in the 
variable z. The generalized method of Madsen et al. (2002, 2003) allows to find 
exact infinite series solutions instead of approximate solutions. 
We first introduce two arbitrary expansion levels z\ and z^ in each layer, namely 
Z\{X) = —o~ih{X) with < (Ti < a and Zi(X) = —a 2 h(X) with a < er 2 < 1, 
and the associated unknowns <f>i and Wi such that 4>i = 4>i{X,z = Zi(X)) and 
Wi = Wi(X,z = Zi(X)) for i g {1,2}. We now look for solutions to the Laplace 
equations in the form of infinite Taylor series in the vertical variable (z — z,), 

4 (X,z) = ^(z-z 4 )> ( r ) W, (3.2) 

n>0 

where the choice of the vertical variable (z — zi) instead of z actually allows to save 
one step compared to the procedure of Madsen et al. (2002, 2003). Injecting (3.2) 
into the Laplace equations (3.1) and using the mild-slope assumption leads to the 
recurrence relation 

A0 ( f 5 - 2(n + l)Vzj • V0 ( ; l+1) + (n + 2)(n + l)0 ( ," +2) S . (3.3) 
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Observing that A(Vzi ■ * ) ~ Vzi ■ AVcf)^' for all k yields the expression of 
drf n ' and <^/ 2n+1 ^ in terms of <b; and Wi, 



(k) 



i(2") 



,(2n+l) 



(2n)! * 
m (-1)" 



(-1)' 



(2ra - 1 



A*i2i + (-l) n - 



2ri 



(2n + l)! v ' (2n + l)! 

Plugging these expressions into the ansatz (3.2) provides the desired expressions 
for the velocity potentials <pi and vertical velocity w, in terms of 4>i and it), for all 
(X, z) belonging to fi,, 



J <fo(Af, z) = C(z - Zi)<t>i +S(z - + Vzi ■ , 
| w t (X, z) = —S(z- Zi) + C(z - Zi)Wi + Vzi ■ T Wi , 
where C and S are infinite-series pseudo-differential operators defined by 

C(A) = V(-l) n — -A™ , 5(A) = V(-l) n -^ -A n , 

v ; ^ ' (2n)l ' v ' ^ K ' (2n+l) 

n>0 V ; ri>0 v ; 

and where the slope terms IV and T Wi are given by 



r,/ 
r„ 



(Z - Zi) 

{z - Zi) 



C(z - Zi)W<t>i +S(z - Zi)Vwi - S(z - Zi)S74>i , 
— S(z — Zi)\?Atfii + C(z — Zi)XJwi + S(z — Zi)Vwi 



(3.4) 



(3.5) 



(3.6) 



The expression (3.4) of provides a theoretical formulation of an exact solution to 
the Laplace equations (3.1). Strictly speaking, we can verify that they are in fact 
solutions to (3.1) with residuals of order 0(\Vh\ 2 , Ah) which are negligible within 
our mild-slope approximation framework. 



(b) Truncation of the Taylor series 

The previous solutions to the Laplace equations (3.1) are purely theoretical 
since they involve infinite-series pseudo-differential operators. To deal with this 
problem, we obviously need to truncate the series at a finite order, and this raises 
the question of choosing the order of truncation. Through this choice, we have 
to reach a compromise between the accuracy of the truncated expression and the 
numerical complexity of the final model. In fact, the truncation order is a key factor 
for the domain of validity of the model: the higher the truncation order is, the better 
dispersive effects are reproduced, so that the model is applicable in deeper water. 
We recover here the common paradigm encountered in works based on asymptotic 
expansions of the solutions to (3.1) where retaining higher order terms increases 
the domain of validity in intermediate or deep water. 

Within our double-layer framework, the increase of the number of unknowns 
allows us to lower the truncation order in comparison with the works of Madsen 
et al. (2002, 2003), and even of Jamois et al. (2006). We take advantage of this 
possibility, and truncate the operators C and S by retaining only the first two 
terms of the series, which leads to the following approximations 

C(A) = 1- yA + (9(A 4 A 2 ) , S(A) = A-yA + 0(A 5 A 2 ). (3.7) 
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We first plug these approximations into (3.4), which leads to the following truncated 
expressions of ft an d u>. 



' tfn{X, z) = (1 - OLiA)4>i + (ft - j l A)w l + Vzi ■ ft ((V - a,VA)ft 
+(ft V - 7,VA)w ? ;) - (ft V - 7*VA)ft] , 
Wi(X,z) = (-ftA + 7 l A 2 )ft + (l-a l A)u; l +Vz 4 - ft((-ftVA 



(3.8) 



-T^A 2 )^,; + (V - aiVA)wi) + (ftV - 7i VA)w 



where 



di{z) = ^ Z ^ , f} l (z) = z-z l , ji(z) = 



(3.9) 

We now reformulate these expressions by applying the operators Pi , defined for all 
smooth enough scalar-valued function u by Pj u = u — ft V 1 Zi ■ Vu, which yields the 
following expressions for all {X, z) in fij , 



Pi<j>i ss (1 - OjA)^,- + (ft - 7iA)ii5i - V* ■ (ftV- 7l VA)0, 



PiWi 



(-ftA + 7i A 2 )ft + (1 - aiA)un + Vz, • (ft V - 7;VA)iS 



(3.10) 



The advantage of this new formulation will become clear in the sequel. 

Remark 3.1. The chosen order of truncation can be motivated as in remark 2.2 by 
a dimensional analysis in the case of a flat bottom. Scaling here z and the expansion 
levels Zi with h\ = ahc, in the upper-layer and /12 = (1 — o~)ho in the lower-layer, the 
dimensionless form of (3.4)-(3.5) is obtained by replacing A 2 by fa in (3.5), where 
Hi = hf / X 2 corresponds to the dispersion parameter in each layer. Supposing that all 
the derivatives are of order 0(1), we can analyze the order of the third terms of each 
series in deep water, for instance kho = 10, where fa is considerably higher than in 
shallow water. If we restrict a to the range [0.25,0.75], this leads to the following 
estimates for n = 2 in deep water: ji?/((2n)!) w 0.08, Mi7(( 2 ™ + 1)!) w 0.01 and 
/i™ +1 /((2n+l)!) ~ 0.02. Thus, the third term of each infinite series in (3.4)-(3.5) is 



very small in deep water: these terms and all the subsequent ones can be neglected. 
The same kind of analysis performed for an uneven bottom leads to the same result. 
This motivates the truncation order chosen in (3.8). 



(c) Pade approximants 

Using the truncated expression (3.10), it is possible to construct an approxi- 
mation to the static Dirichlet-Neumann operator, but involving up to fourth-order 
differential operators at first order in h and fifth-order differential operators in the 
slope terms. Consequently, we now present a method to lower the maximum order 
of the derivatives in (3.10) based on Pade approximants. We follow the strategy 
introduced by Madsen et al. (2002, 2003) and expand the variables ft and Wi in 
terms of auxiliary variables ft and w* through the relations 

& = M i (ziV)<j>t, w i = M i (z i V)w* , 
M»(z<V) = 1+PizfA + qi ZiVzi ■ V , 
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where pi and qi are arbitrary coefficients to be determined. 

We now plug this ansatz into (3.10) and conduct the same dimensional analysis as 
previously, which yields the following expressions: 

Pi<k S (1 - (<* - Pi zf)A) # + fa - ( 7i - PiPizf)A) w* 



+Vz t ■ 



((qiZi - Pi)V + (7, - PiPizf - a 4 (^ + 4p z )z t )VA) 



+ (PiQiZiV - ji(qi + 4pi)ziVA) 

J (3-12) 
\ Wl « (-&A + ( 7i - ft m 2 )A 2 ) 0* + (1 - (a, - Pi^)A) < 

+Vz 4 • [ (-pi( qi + iPi)ziVA) (f>* + ((ft + q tZi )V 
-(li - PiPizf + cti(qi + 4p,)z,)VA) w* 

where we have again kept the first two terms in each modified truncated scries, 
except the fifth-derivative of <p* appearing in the slope terms of PiWi. This choice 
is motivated by the mild-slope approximation, but it clearly unbalances the global 
structure of the slope terms and does have a negative impact on the linear shoaling. 
However, we present in §4 a remedy to this problem. 

We now aim at lowering the maximum order of the derivatives in (3.12) while 
preserving the overall accuracy of the truncated expressions (3.10). This goal can 
be achieved by choosing the coefficients pi and qi in order to introduce Pade ap- 
proximants in the equations. In Madsen et al. (2002, 2003), the authors use Pade 
approximants as a way to improve truncation accuracy without increasing the order 
of the derivatives. In the present work, we rather view the Pade approximants as 
a means to cancel high-order derivatives in (3.12) while preserving the accuracy of 
the truncated expressions (3.10). 

Practically, lowering the maximum derivative order in (3.12) means requiring 
that the factors 7, — PiPizf and qi + 4pi (respectively in front of the fourth-order 
and third-order derivatives) vanish in each layer, via an appropriate choice of the 
constants pi, qi, and Ui (that define the expansion levels z{). Since the quantity 
7i — PiPi%i depends on the vertical variable z through 7, and Pi, it is impossible that 
this factor vanishes over the whole still-water column. Nevertheless, the truncated 
expressions (3.10) of Pi<f>i and PiWi need to be evaluated only at some levels, namely 
z = for the upper boundary of f2i, z = z for the interface and z = z for the sea 
bottom. Consequently, requiring that the quantity 71 — P\p\z\ vanishes at the still- 
water level and at the interface, and that the quantity 72 — P2P2Z2 vanishes at the 
interface and at the sea bottom eliminates all the fourth-order derivatives from the 
model. Using (3.9), this yields 



1 fl-G 



Pi = 7 > P2 = 7 



6 ' r 6 \\ + a 
a (7 + 1 



<ri = -, °z 



2 , 2 /l - (T N 2 



Then, taking qi = —Api yields q\ = and q-2 = — ■ We observe that 

3 3 \1 + a J 

the expansion levels z\ and Z2, respectively defined by z\ = —o\h and Z2 = —a^h, 
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are thus taken respectively at the middle of the upper layer at rest and at the 
middle of the lower one. 

We can now plug the previous values for p^, q i7 and z^ into the expressions (3.12) 
of Pi4>i and PiWi, and evaluate them at the three boundary levels z — 0, z = z, and 
z = z. We first obtain the following expressions at the still-water level: 



1 + ■ V ) 0o = (1 - a*A) 0; + Plw\ + Vh • 7l * V0^ - <5J V< 



(3.13) 



where 



«i = ^ 2 , # = 7? = ?U, <S = ?JU 2 , *i = ^- (3-14) 



At the interface z = z, we use the mild-slope approximation to obtain 
(1 - f#VA -V) 0~I « (1 - a* A) <j>{ - /3>* + Vh ■ [-eJV# + ^ Vw*], 
(1 - §/?JV/i -V) S /9f A# + (1 - a* A) + V/i • [7*Vw;], 
(l + £±i^ V ft. V) ^ « (1 - a* 2 A) 0* + f3*w* 2 + Vh- [ 72 *V0* - 5*Vu, 2 *], 

1 + 2±i^V/i -V) 2^ « -/3aA02 + (1 - a2 A )^2 - V/i-^VwI], 
where 

* ( 1 - cr ) 2 ,2 „* i- 0- , * 5o " + 1 n 
a 2 = ^ — h > P2 = —^— h , 72 = 12 (1 ~ <r)h , 

^2 = JTj ^ . £ 2 = ^-(! -0-)ft • 

Finally, the expressions at the sea bottom z = z are 

1 - ^/3* 2 Vh ■ v)^ S (1 - a*A) 0* - + Vh ■ [S* 2 Vw* 2 - £ |V^] ; 
1 _ £±1/£VZ» ■ v)«* « /?2 A02 + (1 - "2 A) + V/i • bSVufl. 



(3.15) 



(3.16) 



(3.17) 



( d) Formulation of the approximate static Dirichlet- Neumann operator 

The final step in deriving our approximate static Dirichlet-Neumann operator 
Go PP [h] is to reformulate the boundary conditions at z = z and z = z in terms of 
Pi4>i and PiWi instead of 0j and Wi. At the interface, the two operators acting on 
0i and wi on one side, and 02 and W2 on the other side, differ from each other. A 
simple reformulation of the continuity condition 0i = 2 is 

1 - IftVh • V) Ti » (l - ~Wi • V) (l + Z±± /3* 2 Vh -V)fc. 
Consequently, the continuity conditions (2.5e), (2.5f) take the form 



1 - |/3J Vh ■ V) 0i = (l - ^Vh ■ v) (l + Z±±fiVh ■ V)0 2 , 
1 - I^V/i- v) 2JI = (l- ^V/i- v) (l + -^y-^Vft- v)u£ 



(3.18) 
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Finally, applying the operator 1 — ((er+l)/2)/?2 V/i- V to (2.5g) leads to the following 
reformulation of the kinematic condition at the bottom 

(i-^^p;vh-v^m + vh-v(i-^-^-(3;vh-v^ S o. (3.19) 



Gathering all the previous results, we are able to construct the system of five 
equations on the six unknowns (f>o, wq, <f>*, id*, (j) 2 , and w 2 that defines the ap- 
proximate operator Qq PP [K\ linking wq to <f>Q. The first one corresponds to the re- 
formulated Dirichlct condition (Pi4>i)\ z =o — {Pi\z=o)4>Oi l - c - the first equation of 
(3.13). The second and third equations correspond to the continuity conditions at 
the interface z = z, and are obtained by plugging (3.15) into (3.18). The fourth 
one is the condition at the sea bottom derived by plugging (3.17) into (3.19), and 
the last one corresponds to the Neumann condition (PiWi)\ z —o = {P\_\z=o)u>o & s 
expressed by the second equation of (3.13). These five equations can be recast as 





f Mn 


M12 









( ^ ^ 




( Pofo 


\ 




Mil 


M' 22 M' 23 


M' 24 






< 











M 3 i 


M' 32 M' 33 


M' 3i 






02 











{ o 


M' 43 


M' u 


) 








I o 


/ 


P W = 


-#A# + (1 


-aJA 











with the differential operators 



' = 1 + -RVh • V, 
Mn = 1 - a* A + 7* Vh • V , Mu = A* - tffVft • V , 

M 2 i = 1 - a* A - e*V/i • V , M' 22 = -PI + 8{Vh ■ V , 

M' 23 = -(1 - a*A) - Vh ■ ((7 2 * - £)V + ^VA) , 

M' 32 = l- a* A + 7*Vft • V , M' 33 = p* 2 A - ~/3 2 *V/i ■ VA , 

M' 3i = -(1 - a*A) + Vh ■ ((e* + ~)V - ^a|VA) , 
, M' i3 = /3 2 *A + V/i • (V - a*VA) , M' u = 1 - a* 2 A + (7* - /3*)V/i ■ V . 



Incidentally, we observe that it is possible to eliminate all the third-order derivatives 
from the operators M' 23 , M' 33 , M' 3i , and M' 43 in these equations. Indeed, applying 
the operator V/i • V to respectively the fourth, first, and second equation of (3.20), 
using the mild-slope approximation and combining the results yields 



a* 2 Vh ■ VAw* S V/i • Vw* 2 + (3* 2 Vh ■ V A(f>* 2 , 
aJV/i • VA0* S V/i • V</>* + ftVh ■ Vwl - V/i • V^ , (3.22) 
a* 2 Vh ■ VA02 S 2/3* V/i • Vw* + V/i • V0*. + /3* V/i ■ Vwj - V/i • V0 O ■ 



Article submitted to Royal Society 



A double-layer Boussinesq-type model 



15 



Plugging these results into (3.20), we formulate our approximate static Dirichlet- 
Neumann operator Qq PP [K\ as follows: 



/ Mix Mu \ / <j>l \ 

M 2 1 M22 M.2Z M 2 4 

M31 M32 M33 M34 

M42 M43 M44 J \ w% J 



\ 

Pqw 



( Po<t>o \ 
Qi4>o 

Q2<t>0 



(3.23) 



-ptAcft + (1 - a\A - e^Wh ■ V) 



where the differential operators Po, Mu, M12, M21, and M31 are given by (3.21) 
and the new operators are defined by 



M22 
M24 

M33 
M42 

Qi = 



-PI + s;wh • v, 

PI A 



V/i- V, 

1 — a 

-2p{Vh ■ V, M43 = P* 2 A, M44 = 
-V/i - V , Q 2 = -Vft-V, Q 3 

4 cr — 1 



■M 2 3 = -(l-c^A)-7*V/i-V, 

X32 = 1 - at A + (7* - -^-h)Vh ■ V, 

1 — a 

M34 = -(1 - a* 2 A) + (e* - h)Vh • V, 



1 - a* 2 A + (7* - 2p*)Vh ■ V, 
= -V/i- V. 



(3.24) 

We denote by .M = (Mij)i<ij<4 the matrix differential operator linking 0*, 
•j*, 02, and w 2 to 0o in (3.23). We can write each of the operators Mij as a sum 



of a first-order (in h) operator Vij and a mild-slope operator Vh ■ Qij , yielding 



M =P + Vh- Q 



(3.25) 



where 



V 



\ 



Q = 



piA 



' 7i*V 
-eJV 






A 


01 







° \ 


A 


-A* 


-l + o^A 




-PI 




1 - a<*A 


/3*A 


-1 


+ a* 2 A 







P*A 


1 


-a* 2 A J 











^ 


(St 


- |/3i*)V 


-72V 




5* 2 V 


(7i* 




1— (7 


(£2 






-2/3*V 





(72* 


- 2/? 2 *)v y 



We denote by i7 the vector (<f>\, w\* , 4> 2 ,w 2 ) T an( i by ^ the right-hand side 
(Po0Oj Qi0Oi Q2<t>0i Q30o) T - The differential system in (3.23) then takes the form 

OP + V/i- Q){/ = F. (3.26) 

A Fourier analysis of the differential operator V shows that it can be inverted in the 
case of a flat bottom. Considering an uneven bottom, we can write V = Vq + Vh 
where Vq = V(h ), h being the mean depth, and where Vr is of order 0(\7h). This 
ensures that the differential operator V is invcrtible for V/i small enough. Finally, 
the differential operator 1Z defined by 



K = T- 1 {Id-Vh- QP- 1 ) , 



(3.27) 
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where Id is the identity operator, is thus an approximate inverse of the operator 
V + V/i • Q up to 0{\Vh\ 2 ) terms. Hence, 



U « 1ZF 



(3.28) 



which yields the explicit expressions of 4>l , id* , 0*, , and in terms of i 
The very last step consists in introducing the operators 



'o- 



M = -A* A , 

A/2 = 1 - a*A - e*V/i ■ V , 
Qo = 1 - ~#Vfc • V « p - ] 



(3.29) 



and plugging the expressions of <j)\ and w* into the last equation of (3.23), so as to 
obtain the explicit relation between wq, and </>o, and thus the explicit expression of 
our approximate static Dirichlct-Ncumann operator Qq PP [h] 



QT [h] = Qc 



/ M 












( P ° 







N 2 








K 


Qi 
















Qi 




V o 








o J 






) 



(3.30) 



This expression completes the formulation (2.10) of our double-layer Boussinesq- 
type model and will be further improved in §4/ to tighten the model shoaling 
properties. Once again, we stress that the major advantage of the operator Qq PP [h] 
is that it is static. Hence, wc can construct it at t — once and for all. 

Once we have computed 0o, wc can compute the values for </> 1; 4>2, and w| using 
(3.28). Therefore, we can recover the vertical profiles of the velocity potentials <j)\, <f>2 
and the vertical velocities w\ , W2 over the whole water column using a generalization 
of (2.9) for any z € (0, rj), and plugging the computed values of Pi and qi into (3.12). 
We point out that we have neglected the third- and fourth-order derivatives in the 
following expressions, to obtain only second-order derivatives. Of course, there is a 
price to pay for this choice, as discussed in the linear analysis of the vertical profiles 
in §4e. We use the expressions 



2 3 

<M*, X, z) = (1 - yA) 0o + (z - yA) w , 

z 2 

w 1 (t,X,z) = -zA^o + (1 — —A) w , 
in the vertical region z € (0, 77) and the expressions 



'Mt,*, z) » (1 + ftV*i • V) [(l - aj(z)A) 0* + (/?*(*) - 7 *(z)a) 

-/3f(z)Vz 4 -V^*] , 
Wi(t,X,z) S (1 + AVzj • V) [-/3|(z)A0* + (l - aJ(z)A) < 

+/3|(z)Vz l -V<l , 



(3.31) 



(3.32) 
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for z G [z, min(0, rj)] if i = 1 and for z G [z, z] if i = 2, where 

z cr 2 1 (T — ct) 2 

«!(*) = 2^ + ^ + J2 h2 ' = 2 (Z + + ^ + 12 2 ' 

< p{{z) = z+ U -h, pt {z)=z + Z+± h> (3.33) 

k 7*0) = f + f *0(* + ^) , tIO) = Jo + + ^y^)0 + • 

Using (3.31) in the region between the still-water level and the free surface and 
(3.32) elsewhere instead of applying (3.32) everywhere seems to provide a more 
accurate description of the nonlinear profiles, as specified by Madsen et al. (2002, 
2003). This property has also been observed during the nonlinear simulations per- 
formed on the present model in §5. 



4. Linear analysis of the double-layer model 

The goal of this section is to analyze the linear properties of the model (namely the 
phase and group velocities, the vertical profiles of velocity potential and vertical 
velocity, and the linear shoaling) and to optimize their accuracy in relation to the 
results of Stokes linear theory. 



(a) Linearization of the governing equations 



In order to investigate these linear properties, we restrict the analysis to the one- 
dimensional case. We linearize the governing equations (2.10) around steady-state, 
which yields 4>\ = <fio and uT\ = wo, and leads to the linearized model 



d t 4>o + gv = , 
<9 f ?7 - w = , 



(4.1a) 
(4.1b) 
(4.1c) 



or equivalcntly, up to 0(1 i 2 r ) terms, 

dt<f>o + 9V = , 
dtV - wq = , 

/ M n Ma 

M21 M22 M23 

M31 M 32 M33 

V M i2 M43 

p w Q = -ftditi 





M24 
M34 
M44 





( ft 


\ 




( 












Ok 


bo 




4>2 






Q 2 < 


bo 


) 


\ W* 2 


J 




^ Q3C 


bo 1 



[1 - a\dl - e\h x d x 



(4.2a) 
(4.2b) 

(4.2c) 

(4.2d) 



where h x is the bottom slope and the differential operators Mij are left unchanged, 
except that they are here 1-D operators. We point out that in this linearized model, 
the slope terms arc kept in order to investigate the linear shoaling properties. For 
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convenience, we apply the operator Pq to equations (4.2a), (4.2b) to obtain 



d t N 
( Mn 
M21 
M31 


W = - 



gN = 0, 
W° = , 
M12 
M22 
M32 
M i2 





M23 
M33 
M43 

hd- 





M24 
M44 

atdl - 



\ 


( 


01 ^ 




( 


00 \ 


















<i>2 






Q24P 


/ 




W*2 ) 




\ 


Qz4>° 1 



~-ih x d x 



Wi 



(4.3a) 
(4.3b) 

(4.3c) 

(4.3d) 



where 



-- P cj>o,N = P QV , Wo 
(a/2)Plh x d x into 4> 



PqWq, and where we have plugged the relation 



P 



-1 jO 



to obtain Qi 



QiP -V° * Qi4>" for 



i G {1, 3}, since Qi are differential operators of order 0(h x ). 

Plugging the expression of <p° and W° in terms of </>* and wl (given respectively by 
the first line of the differential system (4.3c) and by (4.3d)) into equations (4.3a) 
and (4.3b) leads to the final reformulation of the linearized model 



Mn d t <t>*i +Mu dtwl + gN = 0, 
d t N + ft - [1 - a\d 2 x - elh x d x ] wl 

[Mn - QiMn] 4>\ + [M22 - QiM 12 ] wl h 
[M 3 i - Q2M11] 4>\ + [M32 - Q2M12} wl + M33 

-Q$Mll4>\ - Q3M12WI + M 4 3</>2 + M44W2 = 



= o, 

M 2 3< 



¥2 ~ 
. 



M24W2 = 
M34W2 = 



(4.4) 



We now look for solutions of the classical form 

r](x, t) = Ae w , 6 = Lut-kx 
01 = -i{Bi + ih x B 2 )e l6 , 
w\ =i(Ci+ ih x C 2 )e l9 , 
05 = -i(Di + ih x D 2 )e ie , 
wl =i(Ei +ih x E 2 )e l9 , 



(4.5) 



where A, Bi, £?2,Ci,C2, Di, D2, Ei, and E2 are slowly spatially- varying functions 
(i.e. of the general form F(vx) with v « 1), i is the wavenumber, and lu the 
wave frequency. The complex conjugate parts of these expressions have been left 
out for brevity. As stated in Madsen et at (2002, 2003), the B 2 ,C 2 ,D 2 , and E 2 
contributions are necessary because of the bottom slope, and since the velocity 
potential variables are not in phase with the free surface at the lowest order in h x , 
but are so at the next order. 



(6) Linear dispersion relation 

To determine the linear properties of the two-layer model, we substitute the 
desired form of solutions (4.5) into the linear formulation (4.4) and collect terms 
at the lowest order in h x . Thus, we obtain a linear system of five homogeneous 
equations in A, Bi, Ci, D\, and E±. This system has non-trivial solutions if its 
determinant vanishes, which yields the following dispersion relation 

f_ _ jj_ _ 1 + a 2 (kh) 2 + a 4 (fc/i) 4 + a 6 (kh) 6 

gh ~ ghk 2 ~ 1 + b 2 (kh) 2 + 6 4 (fc^) 4 + b e {kh) e + b$(kh) s ' ^ ' ' 
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where c is the wave celerity and where the (ai) and (bi) coefficients are given by 

h = 2S+ T2' &4 = 352 + 5 5+ 144 ' & 6 = S 2 (2S+^), b 8 = S\ 

where S = cr(l — <r)/12. This dispersion relation is compared in §4e to the exact 
linear dispersion relation given by Stokes linear theory, namely 

c 2 oj 2 tanh kh , , „, 

(4.7) 



and to its [6, 8] Padc approximation which has the same rational form as (4.6). 

(c) Linear vertical profiles 

Coming back to the previous linear system in A, B\, C\, D±, and E\, we can 
now express each of the unknowns Si, Ci, D\, and E\ in terms of A, which leads 
to tedious expressions that are not given here for brevity. For a flat bottom, the 
expressions of <pi and Wi on the whole water column are given by (3.32), (3.33) since 

Ml X, z) = (l - a\(z)A] $ + (f3\{z) - ryt(z)A) w* , 

(4.8) 



Wi(t,X,z) = -(3l(z)A<j>* + ^l-a t l {z)A S jw* 



Plugging the ansatz (4.5) for a flat bottom (i.e. without the mild-slope contribu- 
tions) into the previous expressions and using the computed values of B±,Ci, D±, 
and Ei in terms of A leads to the expressions of <fii(z), w\{z) in the upper layer and 
4>2{z) , W2(z) in the lower layer, in terms of k, h, lo, A, and a. Finally, we recover 
the linear vertical profiles over the whole water column using 

. f d>i(z) for £ < z < , . , . f W\(z) for z < z < , , . 

[ 02(-z) for z < z < z , [ w 2 (z) for z < z < z . 

The resulting vertical profiles will be compared in §4 e to the theoretical linear 
profiles 4> s {z) and u> s (z) coming from Stokes linear theory, namely 



(4.10) 



Ag coshfc(z + fr) 

s (z) = — — sm(wi - kx) , 

cj cosn kh 

Agk sinhk(z + h) , 

w s (z) = — '■ — — sin(wf — kx) . 

lu cosh kh 

(d) Linear shoaling 

We now aim at determining the linear shoaling gradient 70 of the double-layer 
model defined by 

t-*t- < 411 » 

In order to determine this shoaling gradient, we use the method proposed by Madscn 
et al. (2002, 2003). Coming back to the linear formulation (4.4) together with the 
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ansatz (4.5), we then collect terms at the next order, i.e. terms proportional to 
the first derivatives of all the variables. Doing this leads to a new inhomogeneous 
system of linear equations on the unknowns A x , B 2 , C 2 , D 2 , and E 2 involving the 
first derivatives of k and h (namely k x and h x ) and the first derivatives of Bi,Ci, Di, 
and Ei . Differentiating the previously computed expressions of B\,C\, Di , and E\ 
in terms of A, the derivatives of B±, Ci, D\, and E\ can be expressed only in terms 
of A x , k x , and h x . Then, differentiating the linear dispersion relation (4.6) allows 
to express k x in terms of k, h, h x , and a. Plugging all these relations into the 
inhomogeneous system on A x , B 2 ,C 2 , D 2 , and E 2 , we are able to eliminate all the 
unknowns but A x and express it in terms of A, h x , kh, and a, thereby yielding 
the linear shoaling gradient 70 of our double-layer model. The detailed analytic 
expression for 70 is not reported here for brevity. 

The computed shoaling gradient will be compared in §4 / to the exact shoaling 
gradient 7,, which was derived by Madsen & S0rensen (1992) using energy flux 
conservation combined with Stokes linear theory, namely 

_ 2fchsinh2fc/i + 2fc 2 /i 2 (l-cosh2fc/i) 
7s ~ (2kh + sinh2kh) 2 ' ( ' ' 



(e) Optimization of linear properties 

The goal is now to optimize the linear properties of our double-layer model 
by minimizing the errors between these properties and the exact linear properties 
coming from Stokes linear theory. To this end, we tune the free parameter a which 
defines the interface level z = —ah. 

The different errors between the model linear properties and the theoretical ones 
are computed as follows. We respectively measure the errors on the phase celerity, 
the vertical profiles of the velocity potential and the vertical velocity, and the linear 
shoaling gradient as 




S a (K, a) = \l ^ — kh) d(kh) , (4.13) 

with a G {c, </>, w, 7}, K being a reference relative water depth, and 



f° (w(z)-W s (z)\ 2 2 _ 2 



Ei(a, kh) = -J^ f v ^ (Q) SV ' ) dz , E*(<r, kh) = ( 7o - 7s ) 



(4.14) 



where c, c s , (<f),w), (4> s ,w s ), 70, and 7 S come respectively from (4.6), (4.7), (4.9), 
(4.10), (4.11), and (4.12). 

We point out that in all these errors, the weighting by -gr t helps keeping the errors 
to a minimum for low wave numbers, like in Madsen et al. (2002, 2003). Doing this, 
we sacrifice some accuracy at very high wavenumbers (i.e. in deep water), but we 
reinforce the model accuracy in shallow water. This weighting by -X- is especially 
well-suited for the shoaling gradient errors, which are far more critical in shallow 
water than in deep water. 
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Table 1. Optimal values for a 



K 


tt/2 


7T 


2tt 


10 


0~opt 


0.473 


0.428 


0.365 


0.314 



At this point, we could have minimized each of these errors individually, but 
doing this leads to quite different optimal values for a, ranging from 0.2 to 0.5. 
Furthermore, the minimization of the shoaling gradient error is quite problematic: 
we will see later that whatever value we choose for a, the range of validity in kh 
is limited. We thus choose to minimize £ c ,£ ( j )1 and £ w simultaneously to infer the 
optimal value for c, and then optimize the shoaling gradient error £ 7 differently 
We start with the minimization of the errors £ c ,£^, and £ w through the average 
error £ to t (K, a) = ^(£ c (K,o~) +£ < / 1 (K,o~) +£ w (K,o~)), and we compute the optimal 
value of a for several typical values of K: the shallow-water value K = tt/2, the 
intermediate depth value K = tt, and the deep water values K = 2tt and K = 10. 
In this work, we do not optimize a for larger values as the vertical profiles have 
systematically shown an error peak of at least 2% within the range kh G [0, 10] for 
larger K, for instance K = 15 or K = 20. Table 1 summarizes the optimal values 
a pt for each value of K. 

Figure 2 plots c/c s to assess the dispersion error on the phase celerity The 
upper panel compares the errors obtained for each value of 0~ O pt- The price to pay 
for extending the linear range of validity towards deep water values is the growth of 
a small error peak around kh = 3. Indeed, we can see that a 2% error is reached at 
the very deep water value kh = 24 for a op t = 0.473 with a very small 0.01% error 
peak at kh 3, whereas the same error is reached at kh ~ 28 for a op t = 0.314, but 
with a 0.04% error peak at kh w 3. However, such an error is not significant, and the 
overall accuracy of the double-layer model for a op t = 0.314 appears to be excellent 
up to very deep water. In the same way, the lower panel of figure 2 compares the 
error on the phase celerity of our double-layer model with a opt = 0.314 with the 
errors obtained for the Pade [6, 8] approximation, the model of Jamois et al. (2006), 
and the one of Madsen et al. (2002, 2003). We remark that our double-layer model 
accuracy is far better in deep water than what is achieved with the Pade [6, 8] 
approximation and the model of Jamois et al. (2006): a 2% error is reached at the 
very deep water value kh w 28 for the double-layer model, whereas the same error 
is already reached at kh ss 18 for the Pade [6, 8] approximation and kh w 12 for 
the model of Jamois et al. (2006). In comparison with these results, Lynett & Liu 
(2004a, b) showed that their double-layer model reaches the same 2% error (not 
plotted here) at kh = 8. As for the model derived by Madsen et al. (2002, 2003), a 
2% error is reached at kh = 30, i.e. at a slightly greater value than the one reached 
by our model with a — 0.314. 

Figure 3 plots the depth-averaged errors (upper panel) and E w (lower panel) 
on the vertical profiles of the velocity potential and the vertical velocity. We remark 
that the difference between the errors obtained with each value of a opt remains very 
small in shallow water. On the contrary, the benefit for taking a = 0.314 clearly 
appears for both the velocity potential and the vertical velocity in deep water: for 
the vertical velocity profile, a 1% error is reached at kh w 4 and a 2% error at 
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kh 



Figure 2. Comparison of linear phase celerity with the exact Stokes result. The top figures 
compare the errors obtained for our model with a = 0.314 (solid line), a — 0.365 (dashed 
line), a = 0.428 (dash-dotted line), and a — 0.473 (dotted line). The top right figure is 
a zoom on the region kh £ [0, 12]. The bottom figure compares the errors for different 
models: the solid line represents our double-layer model with a = 0.314, the dotted line 
the Pade [6,8] approximation, the dashed line the model of Jamois et al. (2006), and the 
dash-dotted line the model of Madsen et al. (2002, 2003). 



kh w 8. As far as the velocity potential is concerned, a 1% error is reached at 
k w 6.5 and a 2% error is reached at kh ss 10. By comparison, the model derived 
by Madsen et al. (2002, 2003) yields a 2% error at kh = 12 for both horizontal 
and vertical velocity profiles. The difference between the errors on the velocity 
potential and the vertical velocity component can be attributed to the fact that we 
have neglected the fourth-order derivative term in the expression (4.8) of Wi(z). As 
mentioned earlier, this choice entails to sacrifice some accuracy on the profile of the 
vertical velocity. Nevertheless, the global accuracy for both vertical profiles is still 
very good, up to the deep water value kh = 10. 

This analysis of the phase celerity and velocity profile errors does not exhibit any 
major advantage for the choice a = 0.473 instead of a — 0.314 in shallow water. 
The additional errors made with a = 0.314 are at most of order 0.2% in shallow 
and intermediate water. On the other hand, the advantage of this value clearly 
appears in deep water as it appreciably extends the linear range of validity of the 
model, especially for the vertical profiles. Therefore, we decide to adopt the value 
(7 = 0.314 in the sequel. 
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Figure 3. Depth-averaged errors on the vertical profile of the velocity potential (top) and 
vertical velocity (bottom) for a — 0.314 (solid line), a — 0.365 (dashed line), a = 0.428 
(dash-dotted line), and a = 0.473 (dotted line). 

(/) Improved model with tightened shoaling properties 

We now consider the linear optimization of the shoaling gradient properties. As 
shown on figure 4, the model linear shoaling gradient 70 (dashed line) only matches 
the exact linear shoaling gradient "f s (solid line) up to kh = 3 and swiftly departs 
from it beyond that value. We have attempted to optimize this shoaling gradient 
individually, but no value of a makes the two curves fit beyond kh = 3. 

Therefore, a more subtle optimization is needed. This can be achieved using 
the following method. Going back to the full formulation of the approximate static 
Dirichlet-Neumann operator Q^ pp [h] (3.23), we introduce a new constant parameter 
r and apply the differential operator 1 + rh\7h ■ V to the last equation of (3.23), 
which yields 



PqWo 



Pi A - rhftVh ■ VA 



1 - a* A - {e{ - rh)Vh ■ V - rha\Vh ■ VA w{ 



(4.15) 



where P n * 



-[3{ +rh)Vh-V 
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This new formulation does not allow further optimization yet, since the additional 
terms cancel out in the shoaling analysis. However, an interesting option is to 
neglect the third-order derivative on w* in the slope terms. By doing this, we 
greatly improve the shoaling gradient properties, as we will see in the sequel. For 
the moment, we rewrite (4.15) as follows: 



PqW = - [3* A - rhPlVh ■ VA 



1 — ot[A. — [e\ — rh)Vh ■ V 



(4.16) 



We then take advantage of (3.22) to eliminate the third-order derivative from the 
previous equation, and obtain 



+ 



6r 
a 

Vh- V 



V/i ■ V 



1 - a* A - (e* + 2rh)Vh ■ V 



(4.17) 



This implies to redefine the approximate operator as follows: 



Ml 












( P o \ 





Af 2 








K 


Qi 








■A/3 





Q2 











) 




{ Q3 J 



(4.18) 



where 



■A/i 
•A/3 



6r 



V/i • V , 7V 2 = 1 



6r 



V/i-V, QS = 1 



- a* A — (ej + 2rh)Vh ■ V , 
r/i) Vi-Va (P *) _1 - 



(4.19) 



Our modified model is thus identical to (2.10) but with Qq PP [h] as redefined above. 
Starting from it, the new linearized model remains essentially the same, except that 
the second equation of (4.4) now is 



d t N* 



Pld 2 x + -h x d x 

a 



1 - a\d 2 x - (el + 2rh)h x d x 



6r 



hxd;, 



0. 



with TV* = PqT]. This new formulation does not modify the phase celerity and the 
vertical profiles, but allows to further minimize the shoaling gradient error. We can 
now compute the new shoaling gradient 70. Optimizing the parameter r so that the 
error £j(K, a) is minimized for K = 10 and a = 0.314 yields r = 0.0076. 

Remark 4.1. A dimensional analysis shows that this small value of r is coherent 
with our earlier choice to neglect the third-order derivatives of w* in (4-16). 

Figure 4 displays the optimized shoaling gradient for the previous value for r. The 
improvement is quite impressive since the new shoaling gradient exhibits a very 
good agreement up to kh w 12. This accuracy is the same as that reached by the 
model of Madsen et al. (2002, 2003). 

To sum up, our final double-layer Boussincsq-type model consists of (2.10) and 
(4.18)-(4.19) with a = 0.314 and r = 0.0076. This model exhibits excellent linear 
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X 




kh 

Figure 4. Shoaling gradient: model with (r = 0.0076; dotted line) and without (r = 0; 
dashed line) optimization with a = 0.314, and exact shoaling gradient (solid line). 



properties: the phase celerity is accurate up to kh = 28, the vertical velocity profiles 
are accurate up to kh = 10 for the velocity potential and up to kh = 8 for the 
vertical velocity component, and the shoaling gradient is accurate up to kh = 12. 
We emphasize that these properties are not affected by a slight variation of a, which 
makes the model robust towards the parameter a. These results are quite similar to 
those obtained by Madscn et al. (2002, 2003), but the main advantage of the present 
model is that it contains lower-order derivatives and fewer equations, especially in 
2DH. The present double-layer approach is hence a very good alternative to the 
most advanced high-order Boussinesq models as it offers almost the same linear 
properties with a lower complexity. 



5. Numerical simulations: nonlinear behaviour 

On the basis of the model (2.10) derived in §2 (we do not use here the version 
derived in §4 / since we consider a flat bottom) , a classical finite-difference scheme 
is developed to study numerically some nonlinear properties of the model in 1DH. 

We consider the propagation of two-dimensional periodic and regular nonlin- 
ear waves, without change of form, over a flat bottom and without any ambient 
flow field. For this situation, numerical reference solutions can be obtained by the 
so-called stream function method, or more precisely the Fourier approximation of 
the stream function (Dean 1965; Riencckcr & Fcnton 1981). Unlike analytical wave 
theories (such as Stokes or cnoidal wave theories), this numerical approach is ap- 
plicable whatever the relative water depth and steepness are, and very accurate 
solutions can be obtained by increasing the number of terms in the Fourier series 
(e.g. 10, 20, or 50 if necessary for very steep waves). This method was previously 
implemented in a software called Stream_HT by one of the authors (Bcnoit et al. 
2002). For the selected application, the domain of interest covers one wave-length 
(L = 64 m ; k = 2-n / L = 0.098 radm -1 ) and periodic conditions are imposed at 
the two lateral boundaries. We consider a still water depth of h = 96 m, so that 
the relative water depth is h/L = 3/2 yielding kh = 3ir m 9.425, which corresponds 
to deep water conditions. The wave height is chosen as H = 6.4 m, so that the 
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Figure 5. Snapshots of the free surface elevation (top) and the free surface velocity potential 
(bottom) at t = 10 T (left) and t = 25 T (right). The solid line is our model for a = 0.314 
and the dotted line is the reference solution computed using the stream function approach. 



steepness is H/L = 0.1 or kH/2 = n/10 as 0.314, i.e. about 70% of the theoretical 
maximum value of the steepness for a stable wave (Williams 1981). These condi- 
tions correspond to highly dispersive and very nonlinear waves. For this case, the 
period computed with the stream function approach (at order 20) is T = 6.094 s, 
yielding a wave celerity of C = L/T as 10.502 ms -1 . The solution obtained with 
the stream function approach for the free surface elevation and the free surface 
potential is imposed as the initial condition in our simulations. 

The time integration scheme is a classical fourth-order four-stage explicit Runge- 
Kutta scheme, which is known to possess a wide stability region. However, owing 
to the nonlinear nature of the considered test case, this scheme can develop some 
high-frequency instabilities. To avoid such instabilities, an 8th-order Savitsky-Golay 
smoothing filter is applied twice after each time step to 77, <f>\ 1 and w\. The price 
to pay for this filter is a negligible loss of accuracy of the model. As far as spatial 
discretization is concerned, all derivative operators in (2.10) are replaced by cen- 
tered fourth-order finite difference approximations combined with periodic bound- 
ary conditions. The covered domain (of one wave-length) is discretized with 32 cells 
of constant size (equal to 2 m) and a time-step of 0.122 s (corresponding to T/50) 
is used during the simulations. We have verified that using a refined mesh of 128 
cells and 200 time steps per period does not yield any significant improvement. 

Numerical integration of the double-layer model (2.10) is performed over a du- 
ration of 25 T. Simulations have been performed with the deep-water (kh = 10) 
optimal value a = 0.314, even if any value in [0.28, 0.36] would lead to very similar 
results. Outside this range, the results quickly deteriorate. Results obtained after 
durations of 10 T and 25 T are plotted on figure 5 for the free surface elevation 
and the free surface velocity potential. The results arc compared to the reference 
solution obtained with the stream function approach (which propagates at constant 
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speed and without change of form). The results appear to be very good since the 
two curves fit very well until t 20 T, after which small discrepancies become ob- 
servable. Since grid convergence has been verified, this difference can be attributed 
to the approximation in (2.9), where we have neglected fourth-order (and higher) 
nonlinear terms. However, we remark that the global aspect of the model curve still 
corresponds to that of the reference solution at t = 25 T. Furthermore, the free 
surface elevation computed with our model only shows a phase shift error with the 
reference solution: the forms are the same and the amplitude of the waves are equal. 
An interesting remark is that we can use these results to compute the nonlinear 
phase celerity error approximatively. Indeed, taking for instance the free surface 
elevation results and measuring the distance between the crests of the two curves 
yields an approximate value of the difference of celerity between these curves. This 
value provides us a measure of the nonlinear celerity error of the model. We found 
that our model with o~ opt = 0.314 exhibits a nonlinear phase celerity error of about 
0.08%, which is an impressive result. To conclude, the model shows an excellent 
nonlinear behaviour, and we can expect its nonlinear range of validity to reach up 
to kh = 10, at least for flat bottom conditions. 
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